####################################################################
# Replication code for
#   Hirose, Imai, and Lyall, "Can civilian attitudes predict insurgent violence?"
#
# Code to produce maps for estimated IRT support and violence 
#
# Author: Kentaro Hirose
####################################################################




########################################################################
library(spdep)
library(maps)         
library(maptools)     
library(sp)           
library(countrycode)
library(PBSmapping)
library(spatstat)
library(foreign)
library(lattice)
library(gdata)
library(raster)
library(endorse)
########################################################################
load("in.village.shp.RData")
load("out.village.shp.RData")
load("provinces.shp.RData")
load("district.shp.RData")
load("cities.shp.RData")
########################################################################




####################################################################################
### Figure 1: spatial distribution of relative support
####################################################################################
pdf(file = "prediction_map.pdf", width = 6, height = 6)
### Logar(5); Kunar(10); Helmand/Hilmand(23); Uruzgan(26); Khost(32)
####################
q <- c(-2, -1.6, -1.2, -0.8, -0.4, 0, 0.4)
COL <- c(rgb(0, 0.2, 1), rgb(0.3, 0.5, 1), rgb(0.5, 0.8, 1), rgb(0.8, 0.9, 1),
         rgb(1, 0.9, 0.8), rgb(1, 0.8, 0.5), rgb(1, 0.5, 0.3), rgb(1, 0.2, 0))
PCH <- 20
PCH.vio <- 20
CEX <- 0.8
cex.vio <- 0.6
cex.vio.all <- 0.4
col.vio <- "black"
border.col = "gray"
border.lwd <- 1
cex.tex <- 0.8
####################
par(mar = c(0, 0, 1, 0))
####################

par(fig=c(0 + 0.5*0.5/2, 0.5 - 0.5*0.5/2, 2/3 - 0.1 + 0.5*(1/3)/2, 1 - 0.1 - 0.5*(1/3)/2))

A <- (provinces.shp$prov_34_na == "Logar")
B <- (in_village.shp$prov10 == 0 & in_village.shp$prov23 == 0 & in_village.shp$prov26 == 0 & in_village.shp$prov32 == 0)
in_sample.shp <- in_village.shp[B, ]
in_sample.shp <- in_sample.shp[, "diff"]

C <- (district.shp$prov_34_na == "Logar")
plot(district.shp[C,], col = rgb(.8, .8, .8, alpha = .05), border = border.col, lwd = border.lwd)

support <- in_sample.shp$diff
if (TRUE %in% (support <= q[1])){
  plot(in_sample.shp[in_sample.shp$diff <= q[1], ], col = COL[1], pch = PCH, cex = CEX, add = T)
}
for (k in 2:7){
  if (TRUE %in% (support > q[k-1] & support <= q[k])){
    plot(in_sample.shp[in_sample.shp$diff > q[k-1] & in_sample.shp$diff <= q[k], ], col = COL[k], pch = PCH, cex = CEX, add = T)
  }
}
if (TRUE %in% (support > q[7])){
  plot(in_sample.shp[in_sample.shp$diff > q[7], ], col = COL[8], pch = PCH, cex = CEX, add = T)
}

title(main = "Logar", font.main = 1)
####################
par(fig=c(0 + 0.5*0.5/2, 0.5 - 0.5*0.5/2, 1/3 + 0.5*(1/3)/2, 2/3 - 0.5*(1/3)/2), new = TRUE)


A <- (provinces.shp$prov_34_na == "Khost")
B <- (in_village.shp$prov10 == 0 & in_village.shp$prov23 == 0 & in_village.shp$prov26 == 0 & in_village.shp$prov32 == 1)
in_sample.shp <- in_village.shp[B, ]
in_sample.shp <- in_sample.shp[, "diff"]

C <- (district.shp$prov_34_na == "Khost")
plot(district.shp[C,], col = rgb(.8, .8, .8, alpha = .05), border = border.col, lwd = border.lwd)

support <- in_sample.shp$diff
if (TRUE %in% (support <= q[1])){
  plot(in_sample.shp[in_sample.shp$diff <= q[1], ], col = COL[1], pch = PCH, cex = CEX, add = T)
}
for (k in 2:7){
  if (TRUE %in% (support > q[k-1] & support <= q[k])){
    plot(in_sample.shp[in_sample.shp$diff > q[k-1] & in_sample.shp$diff <= q[k], ], col = COL[k], pch = PCH, cex = CEX, add = T)
  }
}
if (TRUE %in% (support > q[7])){
  plot(in_sample.shp[in_sample.shp$diff > q[7], ], col = COL[8], pch = PCH, cex = CEX, add = T)
}


title(main = "Khost", font.main = 1)
####################
par(fig=c(0 + 0.4*0.5/2, 0.5 - 0.4*0.5/2, 0 + 0.1 + 0.4*(1/3)/2, 1/3 + 0.1 - 0.4*(1/3)/2), new = TRUE)

A <- (provinces.shp$prov_34_na == "Kunar")
B <- (in_village.shp$prov10 == 1 & in_village.shp$prov23 == 0 & in_village.shp$prov26 == 0 & in_village.shp$prov32 == 0)
in_sample.shp <- in_village.shp[B, ]
in_sample.shp <- in_sample.shp[, "diff"]

C <- (district.shp$prov_34_na == "Kunar")
plot(district.shp[C,], col = rgb(.8, .8, .8, alpha = .05), border = border.col, lwd = border.lwd)


support <- in_sample.shp$diff
if (TRUE %in% (support <= q[1])){
  plot(in_sample.shp[in_sample.shp$diff <= q[1], ], col = COL[1], pch = PCH, cex = CEX, add = T)
}
for (k in 2:7){
  if (TRUE %in% (support > q[k-1] & support <= q[k])){
    plot(in_sample.shp[in_sample.shp$diff > q[k-1] & in_sample.shp$diff <= q[k], ], col = COL[k], pch = PCH, cex = CEX, add = T)
  }
}
if (TRUE %in% (support > q[7])){
  plot(in_sample.shp[in_sample.shp$diff > q[7], ], col = COL[8], pch = PCH, cex = CEX, add = T)
}

title(main = "Kunar", font.main = 1)
####################
par(fig=c(0.35 + 0.4*0.5/2, 0.85 - 0.4*0.5/2, 2/3 + 0.4*(1/3)/2, 1 - 0.4*(1/3)/2), new = TRUE)


A <- (provinces.shp$prov_34_na == "Uruzgan")
B <- (in_village.shp$prov10 == 0 & in_village.shp$prov23 == 0 & in_village.shp$prov26 == 1 & in_village.shp$prov32 == 0)
in_sample.shp <- in_village.shp[B, ]
in_sample.shp <- in_sample.shp[, "diff"]

C <- (district.shp$prov_34_na == "Uruzgan")
plot(district.shp[C,], col = rgb(.8, .8, .8, alpha = .05), border = border.col, lwd = border.lwd)

support <- in_sample.shp$diff
if (TRUE %in% (support <= q[1])){
  plot(in_sample.shp[in_sample.shp$diff <= q[1], ], col = COL[1], pch = PCH, cex = CEX, add = T)
}
for (k in 2:7){
  if (TRUE %in% (support > q[k-1] & support <= q[k])){
    plot(in_sample.shp[in_sample.shp$diff > q[k-1] & in_sample.shp$diff <= q[k], ], col = COL[k], pch = PCH, cex = CEX, add = T)
  }
}
if (TRUE %in% (support > q[7])){
  plot(in_sample.shp[in_sample.shp$diff > q[7], ], col = COL[8], pch = PCH, cex = CEX, add = T)
}


title(main = "Uruzgan", font.main = 1)
####################
par(fig=c(0.35, 0.85, 0, 2/3), new = TRUE)

A <- (provinces.shp$prov_34_na == "Hilmand")
B <- (in_village.shp$prov10 == 0 & in_village.shp$prov23 == 1 & in_village.shp$prov26 == 0 & in_village.shp$prov32 == 0)
in_sample.shp <- in_village.shp[B, ]
in_sample.shp <- in_sample.shp[, "diff"]

C <- (district.shp$prov_34_na == "Hilmand")
plot(district.shp[C,], col = rgb(.8, .8, .8, alpha = .05), border = border.col, lwd = border.lwd)

support <- in_sample.shp$diff
if (TRUE %in% (support <= q[1])){
  plot(in_sample.shp[in_sample.shp$diff <= q[1], ], col = COL[1], pch = PCH, cex = CEX, add = T)
}
for (k in 2:7){
  if (TRUE %in% (support > q[k-1] & support <= q[k])){
    plot(in_sample.shp[in_sample.shp$diff > q[k-1] & in_sample.shp$diff <= q[k], ], col = COL[k], pch = PCH, cex = CEX, add = T)
  }
}
if (TRUE %in% (support > q[7])){
  plot(in_sample.shp[in_sample.shp$diff > q[7], ], col = COL[8], pch = PCH, cex = CEX, add = T)
}

title(main = "Helmand", font.main = 1)
####################
par(fig=c(0.8, 1, 0, 2/3), new = TRUE)

C <- (district.shp$prov_34_na == "Hilmand")
plot(district.shp[C,], col = "white", border = "white")

COL.inv <- COL[8:1] 
legend("right",  
       c("0.4 <", "[0, 0.4]", "[-0.4, 0]", "[-0.8, -0.4]",
         "[-1.2, -0.8]", "[-1.6, -1.2]", "[-2.0, -1.6]", "< -2.0"),
       pch = rep(20,6), col = COL.inv, cex = cex.tex,
       border = NA, bty = "n", pt.cex = CEX, 
       title = "Relative level of \n ISAF support"
       )


lines(c(63.5, 64.5), c(28.5, 28.5), lwd = 1) 
text(64.08, 28.15, "20", cex = cex.tex)
lines(c(63.5, 63.5), c(28.4, 28.6), lwd = 1)
lines(c(64, 64), c(28.4, 28.5), lwd = 1)
text(64.5, 28.85, "40", cex = cex.tex)
lines(c(64.5, 64.5), c(28.5, 28.6), lwd = 1)
text(65, 28.5, "(km)", cex = cex.tex)


dev.off()
#####################################################################################





############################################################
### Figure 6: Location of villages and attacks 
############################################################
pdf(file = "prediction_map_all.pdf", width = 8, height = 8)
############################################################
q <- c(-2, -1.6, -1.2, -0.8, -0.4, 0, 0.4)
COL <- c(rgb(0, 0.2, 1), rgb(0.3, 0.5, 1), rgb(0.5, 0.8, 1), rgb(0.8, 0.9, 1),
         rgb(1, 0.9, 0.8), rgb(1, 0.8, 0.5), rgb(1, 0.5, 0.3), rgb(1, 0.2, 0))
PCH <- 20
PCH.vio <- 20
CEX <- 1
cex.vio <- 0.6
cex.vio.all <- 0.4
col.vio <- "black"
border.col = "darkgray"
border.lwd <- 2
CEX.TEX <- 1
IN.CEX <- 0.475
####################
par(mar = c(0, 0, 1, 0))
####################

plot(provinces.shp, lwd = 1, border = "white")
plot(provinces.shp, lwd = 0.5, col = "white", border = "gray", lwd = 0.5, add = T)

D <- (provinces.shp$prov_34_na == "Logar" | provinces.shp$prov_34_na == "Kunar" | provinces.shp$prov_34_na == "Hilmand" | provinces.shp$prov_34_na == "Uruzgan" | provinces.shp$prov_34_na == "Khost" | provinces.shp$prov_34_na == "Ghazni"| provinces.shp$prov_34_na == "Kandahar"| provinces.shp$prov_34_na == "Laghman"| provinces.shp$prov_34_na == "Nangarhar"| provinces.shp$prov_34_na == "Paktika"| provinces.shp$prov_34_na == "Paktya"| provinces.shp$prov_34_na == "Maydan Wardak"| provinces.shp$prov_34_na == "Zabul")

plot(provinces.shp[D, ], lwd = 0.5, col = rgb(.8, .8, .8, alpha = .15), border = "gray", lwd = 1, add = T)

DD <- (provinces.shp$prov_34_na == "Logar" | provinces.shp$prov_34_na == "Kunar" | provinces.shp$prov_34_na == "Hilmand" | provinces.shp$prov_34_na == "Uruzgan" | provinces.shp$prov_34_na == "Khost")

plot(provinces.shp[DD, ], lwd = 0.5, col = rgb(.8, .8, .8, alpha = .15), border = "black", lwd = 1, lty = 1, add = T)


plot(out_village.shp, col = rgb(1, 0, 0, alpha = .35), pch = PCH.vio, cex = 0.015, add = T)
plot(in_village.shp, col = rgb(0, 0, 1, alpha = .35), pch = PCH.vio, cex = IN.CEX, add = T)


legend(69, 31.7, legend = "In-sample villages", pch = PCH.vio, col = "blue", bty = "n", , cex = CEX.TEX)
legend(69, 31.2, legend = "Out-of-sample villages", pch = PCH.vio, col = "red", bty = "n", , cex = CEX.TEX)


text(73.5, 34, "Logar", cex = CEX.TEX)
lines(c(69.5, 72.5), c(34, 34), lwd = 0.5)

text(73.5, 35, "Kunar", cex = CEX.TEX)
lines(c(71.5, 72.7), c(35, 35), lwd = 0.5)

text(63.5, 28.8, "Helmand", cex = CEX.TEX)
lines(c(63.5, 63.5), c(29.8, 29), lwd = 0.5)

text(66, 38, "Uruzgan", cex = CEX.TEX)
lines(c(66, 66), c(33, 33 + 4.8), lwd = 0.5)

text(72, 33.5, "Khost", cex = CEX.TEX)
lines(c(70, 71.2), c(33.5, 33.5), lwd = 0.5)

plot(cities.shp, add = T, col = "black", pch = 4, cex = 0.7)
cities.name <- c("(Kabul)")
cities.latitude <- c(34.32)
cities.longitude <- c(69.10)
lines(c(cities.longitude, cities.longitude), c(cities.latitude, cities.latitude + 4), lwd = .5)
text(cities.longitude, cities.latitude + 4.2, cities.name, col = "black", cex = CEX.TEX)

lines(c(71, 73), c(29,29), lwd = 2) 
text(72, 28.5, "100", cex = CEX.TEX)
text(73, 29.5, "200", cex = CEX.TEX)
text(74, 29, "(km)", cex = CEX.TEX)
lines(c(71, 71), c(28.8, 29.2), lwd = 2)
lines(c(72, 72), c(28.8, 29), lwd = 2)
lines(c(73, 73), c(29, 29.2), lwd = 2)


dev.off()
#####################################################################################



